Dynamical Matrices and Interatomic-Force Constants from 
Wave- Commensurate Supercells 



Hadley M. Lawler 
Department of Physics, University of Maryland 
College Park, Maryland 20742-4111 

Eric K. Chang 

University of Modena, Via Universita 4, 4^00 Modena, Italy 

Eric L. Shirley 

Optical Technology Division National Institute of Standards and Technology 
100 Bureau Dr., Stop 844 1, Gaithersburg, MD 20899-8441 



Abstract 

We apply standard, first-principles calculations to a complete treatment of lattice dynamics 
in the harmonic approximation. The algorithm makes use of the straightforward "frozen- 
phonon" approach to the calculation of vibrational spectra and addresses some limitations 
of the method. Our prescription's validity is independent of crystal structure. It treats 
polar crystals in a general way, and it incorporates interatomic-force constants reaching 
beyond the supercell that is considered. For a range of materials, phonon dispersion exhibits 
the close agreement with experiment that is now characteristic of first-principles schemes. 
Results for graphite. Si and GaAs are presented. 
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I. Introduction 



Several years ago, density-functional theory (DFT) [1] opened the way to detailed knowl- 
edge of phonon spectra from first principles. The two most common approaches today are 
density- functional perturbation theory (DFPT) [2], where linear corrections to the equilib- 
rium electronic wave functions are calculated with respect to some perturbation, in this 
case ionic displacements, and the "frozen-phonon" method, which solves for the ground 
state with the displacements "frozen" in [3]. The latter is taken up here. 

Direct-method calculations, in conjunction with the Hellman-Feynman force theorem [4], 
have evaluated restoring forces and phonon frequencies by employing two different types 
of supercells. Initial work utilized wave-commensurate supercells, which map the phonon's 
periodicity to the supercell. This allowed for the evaluation of dynamical matrices for a 
very small set of wave vectors where the phonon polarizations were known a priori [5]. 
Later, supercells were designed to compute interplanar forces. Transforming the planar 
forces into the real-space interatomic-forcc constants (IFCs) made possible the expansion 
of the dynamical matrix throughout the Brillouin zone (BZ) [6]. Our strategy returns to 
wave-commensurate supercells, but displaces ions along Cartesian rather than along normal 
coordinates. 

In this respect, our method resembles that of DFPT, where the dynamical matrices 
are derived directly in momentum space, and then interpolated [7]. Accordingly, similar 
techniques can be exploited to account for the longitudinal-transverse (LO-TO) optical 
phonon splitting, which are somewhat different than those often combined with frozen- 
phonon calculations [6, 8, 9, 10] . 

The present scheme, at the level of DFT, exploits the periodicity of the supercells 
and includes all the interactions, while requiring no prior knowledge of phonon modes. 
The computation of planar- force constants discounts interactions extending beyond the 
supercell, even if they are numerically represented in a calculation [11, 12]. This may 
suggest the possibility of a full dispersion calculation with smaller supercells than have 
been necessary until now. 

The method described here proceeds in two steps: first, a set of high-symmetry wave 
vectors, {k}, is chosen, and the corresponding reciprocal-space, interatomic-force constant 
matrices {C(k)} are computed; and second, the counterpart real-space matrices are ex- 
tracted through a fit. A fit is performed, rather than the familiar Fourier transform, to 
conveniently accomodate an irregular wave-vector grid while abiding by crystal symmetry. 
This distinction from the interpolation in Ref. [7] reflects the differing requirements of each 
treatment. Within DFPT, the phonon and electron grids must be equivalent and are typ- 
ically regular, whereas in the present treatment there is freedom to select the wave-vector 
grid, provided that the corresponding periodicities are commensurate with carefully chosen 
supercells. 
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If, at the Bravais lattice site indicated by the vector R, the a^^ coordinate of the r*'^ 
ion's position is given by: 

rra(R) = Ra + Sra + Xrai'R), (1) 



where the first two terms on the right sum to the equilibrium position, then the real-space 
IFCs are defined by: 

5x™(R')9x^'f,(R + R' 

where V is the Born-Oppenheimer energy surface. The translational symmetry of the crystal 
requires that the IFCs are independent of R'. 

The real-space and reciprocal-space matrices are related by: 

C(k) = ^^C(R)e'''-^. (3) 

The above expression establishes the general fit condition on the real-space matrices, and 
gives an expression for the reciprocal-space matrices at arbitrary k once the real-space 
matrices are known. 

The dynamical matrix is obtained by dividing the reciprocal-space IFCs by the ionic 
masses: 

Polar crystals, or those in which sub-lattice shifts are accompanied by macroscopic 
electric dipolcs, exhibit splitting of the small-momentum, LO-TO phonon degeneracy. This 
effect requires additional consideration to be incorporated into the above procedure. In 
this work, such effects are incorporated with a method analogous to that used in earlier 
work, where a distinction is made between "covalent" and dipole-dipole contributions to 
the real-space IFC matrices [7, 13]. 



II. Sampling of the Brillouin Zone 

The periodic ionic displacements entering a frozen-phonon calculation are constructed 
with cells larger than the primitive cell, or supercells. In particular, cells are constructed 
such that ionic positions are given by 

ri-a(R) = Ra + Sra + dra COs(k-R). (5) 

Above, the displacement from equilibrium in Eq. 1, Xra^R), is written to represent a stand- 
ing wave of wave vector k. As an illustration, considering constructions four times the 
size of the primitive cell or smaller, the face-centered cubic lattice can be generated by 



3 



supercells corresponding to inequivalent wave vectors along the following Cartesian direc- 
tions: (1,0,0), (1,1,1), (1,1,0), (2,1,0) and (3,3,1). Lower-symmetry wave-vectors can be 
represented with larger supercells. The Hellman-Feynman force is computed as 



drraCR) 



within a DFT calculation [14]. Electronic states are solved using pseudopotential tech- 
niques [15, 16, 17], specially-chosen points to sample the BZ [18], and well-known methods 

for handling large matrices [19]. 

The harmonic approximation implies the forces and IFCs arc related by: 

FraiB!') = -j:n'r'l>^rr'abi^ -^'K'bC^) (6) 
= -EKV'.C'./a.(RV/6COs(k-(R" +R')), 

where R is shifted, and the ionic displacements from equilibrium are substituted from Eq. 5. 
If d^/jj is nonzero for only a single ion r and direction b, and is varied over two successive 
DFT calculations, the following demonstrates that the finite difference in Hellman-Feynman 
forces can be related directly to the dynamical matrix. Starting with 



= -EC^./a.(R')cos(k.(R" +R')), 



^ ^ R' 



(7) 



and using the basic identity 

cos(a + h) = cos(a) cos(6) — sin(a) sin(6) 

and Eq. 3, we obtain 

Re {C^r'am - -2cos(k.R") \ dd^,, + dd^,, ) 

T (n 1 ( gF^a(R") 9F^a(-R") ^ 
Im (C,,,„,(k)) = ^-^^ 1^^^^;^ . 

The lattice vector R is arbitrary, but chosen within the supercell such that the denominator 
is nonzero. In the case that k is along certain high-symmetry directions and is half of a 
reciprocal-lattice vector, k • R is a multiple of tt for all R . At these points the condition 
C(k) = C(— k) is valid, and by Eq. 3, the equality C(k) = ^C(— k)^ follows, demonstrating 
that the reciprocal-space tensor must be purely real at any such point. 
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The practical result of Eq. 8 is that an entire row of the reciprocal-space matrix is 
known from two calculations, each equivalent to a supercell calculation of total energy. 
Typically, knowledge of the reciprocal-space IFCs at two or three wave vectors along each 
high-symmetry direction is sufficient for a good phonon dispersion throughout the BZ. 
Although complications arise from the fitting to a real-space representation, and the inter- 
polated dynamical matrices do depend on real-space cutoffs, at least at the level of DFT, 
the formalism here incoporates interactions between ions that are not local to a supercell. 
While the standard approximations to the electronic ground state are present in the DFT- 
calculated dynamical matrices, through construction of supercells representing standing 
waves mapped to the BZ point considered, Eq. 9 formally accounts for all the IFCs. 

In order to carry out the interpolation of the dynamical matrix, a discrete Fourier trans- 
form on the complete set of vectors within a parallelepiped box may appear tempting. If one 
wishes to respect the symmetry of the crystal, however, this tactic will not work. Generally, 
the set of vectors spanned by such a box will include some vectors within a symmetry-related 
subset, and not others. If the set is inclusive enough, the artificial asymmetry introduced 
may be negligible, because C(R) is expected to diminish reasonably quickly with R, and the 
set can include all members of any subset whose components have significant magnitude. 

Rather than performing a Fourier transform, by overdetermining the real-space IFC 
matrices used in Eq. 3 and then fitting them, the crystal symmetry is preserved here, and 
no nuisance is posed by the choice of wave- vector sampling. The following section details how 
this procedure is implemented. Key features of the formulation are that all symmetries of the 
crystal are respected, the symmetries are exploited to minimize the number of parameters 
to be fit, and frequencies of the acoustic phonons are guaranteed to go to zero in the limit 
of small wave vector. 

The next section assumes that the real-space matrices are short-ranged, and do not 
include the long-range Coulombic, ion-ion interactions. As a result, the treatment, with- 
out being supplemented, can be applied only to nonpolar crystals. A separate section 
subsequently extracts the contributions from the Coulombic interactions and extends the 
formalism to polar crystals. 

III. Fitting Algorithm 

The Born-Oppenheimer energy surface of a lattice is given by 1/({xt-(R)}), where Xt-(R) 
denotes the displacement from equilibrium of the atom located at R + s,- , or at basis vector 
St- indexed by r in the cell specified by the lattice vector R. The matrices, {C,-t-/(R)}, are 
the 3x3 matrices denoting the IFCs. They satisfy 

F({x,(R)}) = Fo + ^ E x,^(S)-C.(R)-Xr'(S + R) + ---- 

S,R,tt' 
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Because of crystal symmetry, there are simple symmetry relationships between interatomic- 
force constants. Let U„ be a symmetry operator of the crystal group (i.e., a 3 x 3 proper 
or improper rotation matrix) indexed by n, which satisfies, 



UnSr = Sa + gn + R-(Un, S^) 



(9) 



(10) 



R' = U„R + R(U„,s^/ 



)-R(U„,s^), 



(11) 



where gn is the glide associated with U„, is the atomic position of the r atom relative 
to a Bravais lattice site and, R(U„, s^-), R', and R are lattice vectors. The basic symmetry 
transformation relating interatomic force constants is 



where ea and are unit-vectors in Cartesian directions labeled by the indices a and b. 

We may call the triplets (R, r, r') and (R!,a,a') symmetry equivalent if they satisfy 
Eqs. 9, 10, and 11, for some symmetry operation of the crystal group, U„. We may divide 
these triplets into triplet sets in which all of the members within the same set are symmetry 
equivalent. After we divide all the triplets into sets, we assign to each triplet an ordered 
pair, where j labels the set to which the triplet belongs, and i labels which member 

in the set the triplet corresponds to. There is then a one-to-one correspondence between 
the triplet (R, r, r') and the ordered pair (j, i). 

We are interested in computing the general form of Ct-t'(R), that is, finding a minimal 
set of 3 X 3 basis matrices, {F^, A; = 1, . . . , M}, in terms of which Crr'(R) can be generally 
expressed as a linear combination. That is, we have 



To find {Ffc, k = 1,. . . , M}, we use the basic symmetry relation Eq. 12, and other elementary 
properties of the matrix Ct-^'(R), to form general constraints on Ct-,-'(R). We express these 
constraints as homogeneous linear equations that give us {F^}. Plainly we have 



{VnBaf ■ C^a'i^') ' (U„e6) = 6^ • C^r'(R) " ^b, 



(12) 



M 




(13) 



k=l 



(14) 



dxT-ih{^)dXra{^) 9x^-0 (0)5x^/6 (R)' 



and 




(15) 



dXra{^)dXr>h{'£i + S) dXra{'i)dXr'b{^) ' 
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Eq. 14 is the equality of mixed partials, and Eq. 15 results from the discrete translational 
symmetry of the crystal. Combining Eq. 14 and Eq. 15 and setting S = — R, we obtain 

Combining the definition in Eq. 2 and Eq. 16, we have 

C^,,(-R) = C,,,(R-)- (17) 

Prom this it follows that Qtt' (k) is a Hermitian matrix, if taken as a ZN x matrix, where 
N is the number of atoms in the unit cell. Because Eq. 17 is equivalent to the Hermiticity 
of the reciprocal-space IFC Ct-t'C^)' ^'"^ refer to Eq. 17 as the "Hermiticity condition." 

For a given triplet (R, r, r'), we define the little group and transpose set (-i/R,r,r' and 
Tr,,t,t') as follows. Lr,t,t' is the subgroup of the crystal group that transforms the triplet 
(R, r, r') into itself. Tr -^^t' is the subset of the crystal group that transforms the triplet 
(R, r, r') to (— R, r',r). T^^T-y may or may not be a group. Prom Eq. 12, we have 

Ct-v(— R) = U^CT-T-'(R)Un,U„ G T[i,t,t'- (18) 

Combining Eq. 18 with the Hermiticity condition (Eq. 17), we have 

C^,,(R) = U^C,,KR)U„,U„ G Tj^,r,r'. (19) 

This constrains the transpose of the real-space IPCs. In particular, if one has cubic symme- 
try, it implies that the diagonal blocks of the dynamical matrix must be symmetric. Further 
constraints can be imposed by the little group: 

Crr' (R) = VlCrr' (R)U, U„ e Ln,r,r' ■ (20) 

Eq. 19 defines the transpose-set constraints, and Eq. 20 defines the little-group constraints. 
The next step is to exploit the above constraints toward finding the basis matrices, {Ffe}. 
Later, further constraints are introduced to enforce the acoustic sum rule. 
We may define another set of 3 x 3 basic matrices, {Bj, i = 1, 2, . . . , 9}, as, 

Bi = XX, B2 = yy, B3 = zz, 
B4 = yz, B5 = xz. Be = xy, 
By = zy, Bs = zx, Bg = yx. 
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Here we employ the dyadic notation whereby zy denotes a 3 x 3 matrix with aU zero 
entries except for a 1 in the 3rd row and 2nd column. We may define the scalar product of 
two 3x3 matrices in terms of the basic matrices, 

{Bi,Bj) = 6^J. (21) 

We define the entries of the 9x9 matrices, S, and M(n), by, 

Sij = {Bi,Bj) (22) 

and 

Mij{n) = {Bi,UlBj\Jn). (23) 
If we express C^t-/(R) as a linear combination of the nine basic matrices, i.e., as 

9 

C^^,(R) =E«^^^' (24) 

i=l 



the little-group constraints (Eq. 20) become 



9 

J2iMij{n) - 6ij)aj = 0, G L^,r,r', (25) 
and the transpose-set constraints (Eq. 19), 

9 

E (Mij (n) - S^,)aj = 0, U„ G Tn,ry ■ (26) 

3=1 

If there are Nl operations in the little group and Nt operations in the transpose set, there 
are a total of A'^' = 9{Nl + Nt) constraint equations, not all of which are independent. Let 
us define a set of A'"' 9-component vectors, {f3^'^\n = 1, . . . ,N'}, such that the system of 
equations, 

9 

E/3f)ai = 0,n = l,...,iV', (27) 

i=l 

is equivalent to the combined system of Eq. 25 and Eq. 26. We then define {v("\n = 
1, . . . ,M'}, with AI' < N' , as a Gramm-Schmidt, orthonormal basis that spans the same 
space as We define a 9 x 9 dimensional projector matrix, P, using 

M' 

P^J = E^'i^^f- (28) 

n=l 



8 



Let {w^''\ k = 1, . . . , M} be the complete set of M orthonormal eigenvectors of P with zero 
eigenvalue. We then have 

F, = ^^.fB,. (29) 

i=l 

In what follows, we consistently use the following symbols to denote the various different 
kinds of indices involved. All of the below indices are integers. 

• /X : indexes the k^-point in the BZ. 

• a,b : are Cartesian vector indices. 

• j : refers to a set of symmetry-equivalent dynamical matrices in real space. 

• i : refers to a member in a set. 

• k : refers to the basis matrices of which the IFCs in real-space are a linear combination. 

• a,a' : refer to the atomic basis within a unit cell. 

• Nh{j) is the number of basis matrices in the j*^ set. 
The variables and quantities which we will be using are: 

• Hji : is the lattice vector corresponding to the j*'* triplet set and i^^ member. 

• "Oi' "^ji • indices of the basis atoms of the j*'* set and i*^ member. 

• C(j) : is a 3 X 3 matrix for the 1** member of the j*^ set. It is the Tji,T'j^ block of 
the IFCs in real-space with separation vector Rji, i.e., C(j) = ^ (R-ji)- 

• Fjfc : for a given j, {F jk, k = 1, . . . , Nh{j)} is a linearly independent set of Nh{j) 3x3 
matrices (the basis matrices). They are related to C(j) by 

Nbif) 

c(i) = E "ifc^i*^' (30) 

k=l 

where the ajk are to be determined by a fitting procedure to be described and are not 
determined by symmetry. 

• C(/x, cr, cr') : is a 3 X 3 matrix giving the a, a' block in the IFCs at Bloch vector k^^. 

• Ujj : is the matrix transformation that relates an element in the set/member ji with 
an element with set/member ji, according to, C7-^..^'.(Rji) = U^C7-^.^^'^(Rji)Uji. 
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We consider the following quantity s, defined as, 

^ = EE I E [Uj:C(i)U,,]„,e*--^- - [C{f,,a,a')U\Ml^), (31) 

aa' abn {i,j\Tji=a,T^^=a'} 

or 

4{ajk}) = E E I E [Uj,F,fcU,i]„ba,fee^'^--^- - [C(m, a, a')U\Ml^), (32) 

aa' abn {i,j,k\Tji=(T,T^^=a'} 

where here we show the explicit dependence on the set of coefficients {ajk} and the weights 
{w{fi)}, which are arbitrary and can be chosen for a variety of convenient ends, such as 
phase-space weighting. The matrix, C(/Lt, o", o"'), is determined by means described in pre- 
ceding sections of this paper. The parameters {ajk} are fitting parameters. Here we choose 
the triplet sets, (j, i), for which we have 

II "^i^ ~^ ji T?^ II ^ 

where Tc is some cutoff selected by the user. 

To ensure that Ct-t'O^) is manifestly Hermitian, we divide the triplet sets into two 
groups. Within group (1) are those for which there exists a paired counterpart in a 
different set, i.e., a pair and such that we have R^j = — Rj'j' and Tji = 

rj,j,,rjj = Tj/j'. Within group (2) arc those for which there exists no such pair. We let 
Q{i, j\Tji = a, rjj = a'} denote the set corresponding to group (2) and P{i, j\Tji = a, rjj = a'} 
denote the set corresponding to group (1) with only one of the pair members counted to 
avoid overcounting. (Note that the triplet sets in group (1) occur in pairs.) This implies 
including half of the triplet sets in group (1) when performing the fit, and using the Hcr- 
miticity condition to deduce Ct-t-/(R) if it belongs to an excluded triplet set. It should not 
matter which of a pair of triplet sets is included in the fit. 

We then rewrite Eq. 32 in terms of these two groups as: 

s{{ajk}) = EEl E ajk[VjiFjkVjiU(^'''-^^'-[C{ii,<7,a%i,\Mf^) 

aa' abn Q{i,j,k\Tji=a,Tr^=(7'} 

+ E E I E a,fc[Uj,F,feU,-,]„5e^^''-^- - [C(/x, a, a')U\Mf^) 

+ EEl E ajk[VjiFjkVjiUe-''^^-^^^-[Ci^i,a\a)U\Mf^). 

aa' abn P{i,j,k\'^ji='^>Tji=<^'} 

We perform a constrained minimization on the quantity s{{ajk}) subject to the acoustic 
sum rule, which can be expressed by 

^[C^^.(R)]ab = 0. (33) 

Rr' 
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Here there is one equation for every combination of r, a, and b, giving a total of 9N 
constraint equations. The acoustic sum rule is therefore a further set of constraint equations 
linear in the fitting parameters. To determine these linear equations explicitly, we first define 
the quantity, 

(34) 



Q{i\Tji=(J,T'.^=a'} 
P{i\Tji=a'y.^=a} 

The acoustic sum rule applies to the wave vector 

k;, = 0. (35) 
It is helpful to take the value of fj, from Eq. 35, and define: 

The acoustic sum rule is satisfied through the fitting parameters which satisfy the following 
system of linear equations: 

ajkCjk-aa'ab = 0. (37) 

jk 

This system of 9N equations is not necessarily linearly independent. For example, in a 
system like diamond, there is only one linearly independent constraint equation in a set of 
18 linear equations. To reduce Eq. 37 to an equivalent smaller set of linearly independent 
equations, we use Gauss-Jordan elimination. Supposing that they reduce to a system of Nc 
linearly independent equations, we have 

5^ {{ajk}) = E =0,iy = l,...,Nc. (38) 

jk 

By minimizing s{{ajk}) subject to the constraint relation in Eq. 38, we have 
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Here A*^''-' are Lagrange multipliers, which need to be found. Eq. 39 reduces to 

jk v=l 

with 

Mj'k',jk= Yl '^it^)^jk;iJ.aa'abA*>k';iM7a'ab (41) 
Hcrcr' ab 

and 

bjk= E <^^)^jk;fUTa'ab[C(.n,a,a')U. (42) 

;ucrcr' 06 

If there are Nf fitting parameters (i.e., the a's), then Eqs. 38 and 40 form a system of 
Nc + Nf equations with Nc + Nf unknowns. 

IV. Macroscopic Fields and Polar Crystals 

For nonpolar crystals, and hence those with no optically active phonons, the method de- 
scribed above is sufficient for the calculation of full phonon dispersion. However, polar 
materials, in which distortions interior to the primitive cell can generate electric dipoles, 
require further consideration. This is because the boundary conditions of DFT calcula- 
tions will not allow any physics associated with macroscopic electric polarizations or fields. 
The latter are connected to the LO modes, and the DFT calculations arc relevant only to 
the acoustic and TO phonon branches. The macroscopic fields can stiffen the small-wave- 
vector, LO modes in a manner that can depend on the direction from which the wave vector 
approaches zero [9]. 

The situation is further complicated by the fact that it is only at the BZ center that 
supercell DFT calculations do not include all the lattice-field interactions. So the fitting 
procedure must reintroduce these interactions at the BZ center while leaving the dynamical 
matrices at other points unaffected. 

To do this, the Born effective charges are first calculated [20]. These charges are defined 

as 

7* — O-^^iL 

Tab — <-i ) 
OUrb 

where P is the macroscopic polarization (in units of dipole per unit volume) and is a 
sublattice shift. These charges must be zero for a nonpolar crystal, and where they are zero 
there is no energy difference between the LO and TO modes in the long-wave limit. 

Considering two ions in the basis r and r' at sites separated by lattice vector R, the 
contribution to the real-space IFCs from the screened electrostatic energy between the two 
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IS 

^-'"''^^^ = dXra{K')dxM^' + ^y ^^^^ 

The Coulomb energy here is associated with dipole-dipole interactions, and is given by: 

1 ^ ^*nh^*'n'h' ( ^aa' o dada' 
n'RTT'bb'aa' 



The above should be understood to exclude terms corresponding to the same ion on the 
same site. The relative coordinate is indicated by d 

d = R + St' — St, 

and 600 is the dielectric constant. The derivatives are with respect to displacements from 
equilibrium along the corresponding ionic coordinates. Gonze and Lee [21] and Shirley et 
al. [22] write in detail the reciprocal-space transform of Eq. 43, 

C^°"'(k) = J2 C^°"'(R)e^''-^ (44) 

R 

in accord with the Ewald formulation [23] . Eq. 44 can be written as: 

C^°"'(k) = C""(k) + C""(k), (45) 

where the analytical and nonanalytical contributions are denoted with superscripts. The 
nonanalytic contribution arises from the G = term in the reciprocal-lattice Ewald sum, 
and accounts entirely for the LO-TO splitting. Its limiting behavior for small wave vector 
can be shown to be [24] : 

r^na /I n^ _ X^aa'^ rabZ ' r' a' b' l^\J>'a' 

- — — — — — . 

The small-k, nonanalytic contribution to C'^°"'(k) is missing in a supercell calculation. 
In order to reintroduce it, Eq. 3 can be modified as: 

C^^^(k) - T(k) = J2 e''''^C(R), (46) 
R 

where C^^'^(k) indicates the matrix calculated within DFT and excludes the effects of the 
macroscopic electric field, and T(k) is defined as: 



^W-i c^o«i(k) : k^^O. 
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Figure 1: Phonon dispersion in Si. 



The distinction between the two cases is always clear, because we are only considering 
the sampled BZ points, and all of these points, except for the T point itself, are far from 
the BZ center. Within the fitting procedure, the real-space and reciprocal-space IFCs are 
now replaced by C(R), and C^^-^(k) — T(k), respectively. 

Then, the full reciprocal-space IFC matrix, with the contribution of the macroscopic 
fields, is: 

C(k) = C^°"'(k) + ^C(R)e^^-^. (48) 

R 

V. Computational Details 

Dispersion curves for Si, GaAs and graphite arc presented in Figs. 1,2,3. Below are 
listed which phonon k points are used in each fit, and which real-space triplet matrices 
are included. These particular selections are somewhat arbitrary, so long as the BZ is 
adequately sampled and the triplet matrices are over deter mined. For the case of graphite, 
extra weight is given to the T point in the fit to ensure that at very small k the very small 
eigenvalues of the dynamical matrix, o;^, do not go negative. 

The symmetries are determined from lattice structure: for GaAs and Si, the lattice 
vectors are 

ai = |(0,l,l), a2 = ^(1,0,1), a3 = ^(1,1,0). 
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Figure 2: Phonon dispersion in GaAs. 
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Figure 3: Phonon dispersion in graphite. 
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A lattice constant of 10.05 ao is employed for Si, and 10.46 ao for GaAs (all reported lattice 
constants correspond to the values for the theoretically relaxed crystals). The basis vectors 
for each crystal are: 

Sri = -g(^i + ^2 + as), = -(ai + a2 + as). 

In Si, the two basis sites are symmetry-mapped, and in GaAs our calculation identifies the 
Ga with Ti. 

The lattice vectors of graphite are: 

ai = a (^-^.^'Oj ,ai = a '^3 = c(0,0,l), 

the in-plane lattice constant is 4.88 ao, and the lattice parameters take the ratio, c/a = 2.72. 
The basis for graphite is 

1 1 1,1 1,1 

Sri = ^as, St-2 = -^as, s^-g = -- (ai -I- &2) + ^aa, 8^-4 = - (ai 32) - -as. 

For GaAs, the k points included in the fit are r, X, L, K, |X, |L, ^X, ^L, fK,W, 

using cubic convention, and the lower-symmetry point ^(f,f,^). The prototypical 
triplets are the basis pairs (tiTi), (rir2), and {T2T2) for R = (0,0,0); the basis pairs 
{tiTi) ■, {T1T2) ■, and (r2r2) for R = |(1,1,0); and the basis pairs (tiTi) and {T2T2) for 
R = a(l,0,0). For Si the DFT calculated phonons are the same set as for GaAs. The 
triplets are pairs (tiTi) and {T1T2) for R = (0, 0, 0); (tiTi) and {T1T2) for R = ai; and (tiTi) 
for R = ai -|- a2 — as. All triplets symmetry-related to those within the prototypical set 
enter the computation. 

The graphite calculation includes the phonons at F, K, M, A, W, |M, ^M, and 
^K. The graphite triplets are: the pairs (tj^ti), (T1T2), (tiTs), {tit^) , (t^t^) , and {T3T4) for 
R = (0, 0, 0); the pairs (t-^ti), (T1T2), (nrs), (tsTs), and (rsTi) for R = -ai; the pairs (nrs) 
and (T2T4) for R = ai — a2; and the pairs (tiTi), (r2r2), (tsTs), and {T4T4) for R = a2 — 2ai. 

VI. Conclusions 

We have made some innovations within the direct approach to the computation of lattice 
dynamics. These include the calculation of dynamical matrices without resort to real-space 
atomic force constants, which may reduce the supercell size necessary for a full dispersion 
calculation; the general extension to polar crystals; and a symmetry-preserving fit algorithm 
for interpolation throughout the BZ. The direct method is shown to render a complete, first- 
principles characterization of the phonon spectrum for a range of materials. Accordingly, the 
direct method can be considered complimentary to density-functional perturbation theory 
as a means of computing a variety of fundamental, phonon-related solid-state properties. 
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